arXiv:1501.00094vl [hep-lat] 31 Dec 2014 


PS 


PROCEEDINGS 

0F SCIENCE 

The CLS 2+1 flavor simulations 



Piotr Korcyl* desy 14-224 

John von Neumann Institute for Computing (NIC), DESY, Platanenallee 6, 15738 Zeuthen, 

Germany 

E-mail: piotr . korcyl@desy . de 

We report on the status of large volume simulations with 2+1 dynamical fermions which are being 
performed by the CLS initiative. The algorithmic details include: open boundary conditions, 
twisted mass reweighting and RHMC, whereas the main feature of the simulation strategy is 
the approach to the physical point along a trajectory of constant trace of the mass matrix. We 
comment on the practical side of the above issues using as examples some of the newly generated 
ensembles, which presently cover lattice spacings between 0.05 fm and 0.11 fm and pion masses 
between 150 MeV and 415 MeV. 


The 32nd International Symposium on Lattice Field Theory 
23-28 June, 2014 

Columbia University New York, NY 


* Speaker. 


(c) Copyright owned by the author(s) under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike Licence. 


http://pos.sissa.it/ 








The CLS 2+1 flavor simulations 


Piotr Korcyl 


In order to suppress or eliminate systematic errors associated with the continuum and chiral 
extrapolations, modern lattice QCD simulations try to explore the parameter space of a ~ 0.05 fm 
and physical pion mass. This is possible not only because of the increase in the available computer 
resources but mainly because of the progress in our understanding of numerical algorithms. New 
algorithmic developments are aimed at alleviating the overall scaling with a and m K and try to 
deal with three phenomena which mainly contribute to the total cost of a simulation. These are: 
the autocorrelations which grow as the continuum limit is approached, partially because of the 
topology freezing; the accidental zero-modes of the Dirac operator which render the simulation 
unstable and finally, the growth of the condition number of the Dirac operator as the light quark 
masses are lowered. The algorithmic improvements trying to circumvent the above issues are the 
following: 

• open boundary conditions in time [2] and periodic boundary conditions in the space direc¬ 
tions, are supposed to let the topological charge flow in and out of the simulated volume and 
thus decrease the autocorrelation times, 

• twisted mass infrared regulator [3] which, by moving the spectrum of the Dirac operator off 
the real axis, decreases the effects of the accidental zero-modes, 

• deflated solver [4] which, makes the cost of the inversion largely independent of the quark 
mass. 

The three mentioned ingredients were implemented in the open-source package openQCD [5] and 
used in large volume simulations with 2+1 flavors of non-perturbatively improved Wilson fermions 
[1]. The simulations were made by a joint effort of the CLS [6] collaboration with computer time 
that was granted through two PRACE projects (at LRZ’s SuperMuc in Munich (Germany) and 
CINECA’s Fermi in Bologna (Italy)) and three national projects (Gauss-Center in Julich (Ger¬ 
many), NIC in Julich (Germany), CSCS in Lugano (Switzerland)). 

In the following we present some of the first impressions gathered during these simulations. 
We start by briefly describing the algorithmic details, then we discuss the strategy to reach the 
physical point, and finally we concentrate on several observables such as to, the topological charge 
evaluated at positive Wilson flow time or the twisted mass and rational approximation reweighting 
factors evaluated on several of our new ensembles (for more details see Re.[l]). 

1. General setup: openQCD-1.2 

The described simulations are being made using the open-source package openQCD version 
1.2 [5]. The package implements the Hybrid Monte Carlo algorithm along with several solvers 
enabling simulations of dynamical fermions. The gauge action used is the tree-level Symanzik 
improved action (co = |, c\ = — ^), while the fermion action is discretized with the Wilson /J(a)- 
improved action with csw determined non-perturbatively [7]. Open boundary conditions for the 
gauge fields are implemented following [8, 2] 

F ot:W| Xo= o = F 0k {x)\ xo=T = 0, k = 1,2,3, 
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for the gauge fields and 

P +VW L=0 = P V'M \x 0 =t = °» P± = ^ l± H>) 

W p -U=o = ^W p +U=r = °> 


for the fermion fields. The twisted mass infrared regulator is introduced for the two degenerate 
light fermions [3] by replacing the original determinant dct(D : D) by i n the simulation 

algorithm and reweighted to the proper distribution at the analysis stage. Furthermore, the later 
determinant is split into typically n ~ 5 — 6 parts and each part is associated to an independent 
pseudofermion field. 


det {D'D + flp 2 t 2 det(ZX>'Z) +ji 0 2 ) det(D 1 D + rf) 

AclUr D + 2 fll) y P " dcl(D'f) + 2p2) det(Z)tD + n* +1 ) 


(1.1) 


The dynamical strange quark is implemented with the RHMC algorithm with reweighting [9]. 
All contributions to the HMC forces are integrated along the trajectory with a set of nested MD 
integrators featuring a 3 level hierarchy (~ 10 steps on the coarsest level, one step on finer levels) 
where 2 nd (level 2) and 4 lh (levels 0 and 1) order Omelyan integrators are used. 


2. Strategy: our way to the physical point 


As in all such projects, one has to choose the starting point in parameter space and define 
the trajectory along which the continuum and chiral limits are taken. To do so, we define two 
dimensionless quantities which are used to perform the matching: 


02 (j 3 , K u , K s ) = 8 t 0 ml oc m ud + 

04(j 3 , Kn Ks) = %t Q (m 2 K + l -m 2 K ) oc tr M + 

Our trajectory is defined by keeping the trace of the mass matrix constant, namely 

^ — = const RMr = const + 0(a) 

Kj 


i=u,d,s 


(2.1) 

( 2 . 2 ) 


(2.3) 


In particular, at the symmetric point K u = K d = k s one is left with only one free parameter and 
therefore the matching for different lattice spacings can be done relatively easily. Figure 1 presents 
the results of the performed matching. The three points on the right (with 02 ~ 0.75) come from the 
three matched ensembles at the symmetric point, whereas the remaining data points have different 
values of 02 but should lie along the line of 0 4 « const. The estimate of the physical point uses the 
values of m K and nix taken from PDG and yTo = 0.1465 fm from [10]. 


3. First impressions 

During the granted computer time 10 ensembles were generated with MC chains of length 
of at least 50 T exp , whose details are summarized in table 1 (for more information see Ref.[l]). 
Simultaneously with the generation of configurations several quantities were monitored, such as 
the topological charge and the action density at positive Wilson flow time and the two reweighting 
factors. We now discuss them in more detail. 
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Figure 1: Tuning to the trajectory of constant trace of the mass matrix. 


m K m K 

0.085 fm 0.064 fm 0.05 fm a [fm] 
3.4 3.55 3.7 [5 

415 MeV 415 MeV 

440 MeV 350 MeV 

470 MeV 280 MeV 

480 MeV 220 MeV 

490 MeV 150 MeV 

32 3 x96 32 3 x96 48 3 x128 

32 3 x96 

32 3 x96 48 3 x128 64 3 x192 

48 3 x96 64 3 x128 

64 3 x128 


Table 1: Overview table of lattice spacings and pion masses. The red entries correspond to the ensembles 
at the symmetric point that were used in the matching between different lattice spacings. The blue entries 
correspond to the ensembles being still generated. 


3.1 Re weighting factors 

Before an average value of any physical observable is computed, one has to measure the 
twisted mass and rational reweigthing factors, used to be denoted by Wo and Wj, respectively. 
They are given by 


_ det(D t D)det(D t D + 2 J u 2 ) 
0 ““ det(DtD + jU 2 ) 2 


Wi = det (D S R), 


(3.1) 


with Z = D]D s R 2 — 1, where det D s = Wj dclR 1 and R is the Zolotarev rational approximation. 
We stress that the inaccuracy of the rational approximation is dealt with in these simulations by 
reweighting instead of an additional Metropolis step as is done typically. 

The left and middle panels of figure 2 show samples of MC histories of the twisted mass 
reweigthing factor for different ensembles. Their values fluctuate close to 1 as expected. From the 
algorithmic point of view there is nothing wrong with having Wo close to or equal to zero, however 
such small reweighting factors may lead to a decrease of the precision of the final average [12], 
therefore they have to be monitored during the entire production phase. 

As far as the rational reweighting factors are concerned, the right panel of figure 2 shows 
extracts of MC histories of W] for the three ensembles at the symmetric point at three different 
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Figure 2: Histories of the reweigthing factors. On the left plot we compare Wo of the three ensembles having 
the same lattice spacing but different m K or /i, the plot in the middle shows the comparison of Wo evaluated 
at the symmetric ensembles at different lattice spacings, whereas the plot on the right shows the histories of 
the rational approximation reweighting factor at the symmetric point. Plotted data comes from ensembles: 
left panel: H105r005, H105r001, C101r014; middle panel: HlOlrOOl, H200r000, N300r002; right panel: 
HlOlrOOl, H200r000, N300r002, see Ref.[l], 



Figure 3: Topological charge at the symmetric point at flow times t/a 2 = 3.0, 6.0, 12.0 which roughly reflect 
the scaling of fo with a. Data comes from ensembles HlOlrOOl, H200r000 and N300r002, see Ref.[l]. 


lattice spacings. We see that the fluctuations are very small and do not depend significantly on the 
lattice spacing. One concludes that the reweighting of the rational approximation works very well. 

3.2 Physical observables 

Once the reweighting factors are computed one can evaluate expectation values of physical ob¬ 
servables. First, we compute the topological charge given by the field-theoretic definition at Wilson 
flow time of the order of fo. The left panel of figure 3 shows the comparison of the history of the 
topological charge at three different lattice spacings evaluated on ensembles at the symmetric point. 
Similar information is plotted on the right panel of that figure, where the integrated autocorrelation 
times are shown. One clearly sees an increase of % mt as a is decreased with a scaling in agreement 
with [2]. 

Next, we discuss the YM action density evaluated at Wilson flow time of the order of fo- Figure 
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•lUIDU l MDU 

Figure 4: The left panel shows the fluctuations of the YM action density at the symmetric point at flow time 
t/a 2 = 3.0, 6.0, 12.0. The right panel shows the comparison of autocorrelation functions of E(t) and Q(t). 
Data comes from ensembles HlOlrOOO, H200r000 and N300r002 (only part of the full statistics is shown), 
see Ref. [ 1 ]. 


4 compares the MC history of the fluctuations of this observable for three different lattice spacings. 
Figure 4 also shows the autocorrelation functions of E{t) corresponding to the two coarsest en¬ 
sembles together with the autocorrelation function of the topological charge. We notice that the 
integrated autocorrelation time of the topological charge is smaller, which is consistent with the 
recent findings of [11]. 

E(t) evaluated at each time-slice independently can be use to exhibit cut-off effects introduced 
by the open boundaries. On figure 5 we plot E(t) for different ensembles as a function of physical 
distance from the boundary expressed in terms of xq jy/to- The value approaches 0.3 by difinition 
in the central region of the lattice. We see that the height of this lattice artifact depends strongly on 
the lattice spacings (the three distinctive sets of data points) and is roughly independent of the pion 
mass (different colors of the data points within each set). We also included results from a simulation 
with a Iwasaki gauge action at a lattice spacing of a « 0.09 fm. The discrepancy between these 
data and the data obtained with the Wilson action with a similar lattice spacing confirms the cut-off 
nature of the visible structure. On the right of figure 5 we show the dependence of E(t,x o) on a 2 
for three different values of xq /y^o- The data cannot be described by a linear dependence on the 
lattice spacing. 

Finally, we present our estimates of to for all the considered ensembles normalized by the 
value obtained at the symmetric point as a function of m K , see figure 6. We also include in the plot 
the analytic prediction from NLO chiral perturbation theory worked out in [13]. We notice that the 
data points deviate from this prediction by very small amounts; cut-off effects are small. 

4. Conclusions 

Several points should be highlighted. The algorithmic improvements turned out to work very 
well in practice. In all our runs the algorithm proved to be stable. The first preliminary measure¬ 
ments of physical observables show that the expected precision can be reached. The trajectory of 
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Figure 5: E(to.xo) as a function of the distance from the boundary shows large cut-off effects. They exhibit 
strong lattice spacing dependence (see plot on the right), but are largely independent of the pion mass. 
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Figure 6: Dependence of to/fl ef on m K at different lattice spacings. The data show that the cut-off effects in 
this observable as well as its pion mass dependence are small, which is favorable in view of the tunning to 
the chiral trajectory. 
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constant trace of the mass matrix allowed precise matching and the tunning to the chiral trajectory 
turned out to be possible with small effort. The autocorrelation time of the topological charge is 
not the largest one, in agreement with Ref.[l 1]. 
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